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Abstract 

The hterature on reconstruction formulas for photoacoustic tomogra- 
phy (PAT) is vast. The various reconstruction formulas differ by used 
measurement devices and geometry on which the data are sampled. In 
standard photoacoustic imaging (PAI), the object under investigation is 
illuminated uniformly. Recently, sectional photoacoustic imaging tech- 
niques, using focusing techniques for initializing and measuring the pres- 
sure along a plane, appeared in the literature. This paper surveys existing 
and provides novel exact reconstruction formulas for sectional photoacous- 
tic imaging. 



1 Introduction 

The literature on reconstruction formulas and back-projection algorithms for 
photoacoustic imaging is vast. Wang et al. developed reconstruction formulas 
for cylindrical, spherical, and planar measurement geometries in a series of pa- 
pers [24, 22, 25], and recently many more algorithms based on reconstruction 
formulas have been developed (see the survey [12]). 

Also different measurement devices for the ultrasound pressure have been 
suggested. Most common are small detectors based on materials, which exhibit 
a strong piezoelectric effect and can be immersed safely in water (i.e. polymers 
such as PVDF). In analytical reconstruction formulas, they are considered point 
detectors. Other experimental setups have been realized with line and area 
detectors, which collect averaged pressure (see [20] for a survey). 



*The work has been supported by the Austrian Science Fund (FWF) within the national 
research network Photoacoustic Imaging in Biology and Medicine, project S10505-N20. 
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Here, we consider the problem of photoacoustic sectional imaging. Opposed 
to standard photoacoustic imaging, where the detectors record sets of two- 
dimensional projection images from which the three-dimensional imaging data 
can be reconstructed, single slice imaging reconstructs a set of two-dimensional 
slices, each by a single scan procedure. The advantages of the latter approach 
are a considerable increase in measurement speed and the possibility to do selec- 
tive plane imaging. In general, this can only be obtained by the cost of decreased 
out-of- plane resolution (i.e. the direction orthogonal to the focusing plane). Ex- 
perimentally, one can obtain photoacoustic sectional imaging by illuminating a 
single plane of the object and by using a focused detector. Technical details are 
provided in Section 2. In our experiments, the measurement data are recorded 
on a cylindrical domain dCl x R, where dfl denotes the boundary of a smooth 
domain f2 in R^. 

The difference in this model to previously studied models is that the wave 
propagation is considered fully three-dimensional, the initialization and mea- 
surements are fully two-dimensional due to the selective plane illumination and 
detection. Therefore, such setups require novel reconstruction formulas. In par- 
ticular, as a further novelty, we present reconstruction formulas in ellipsoidal 
domains. 

This paper surveys existing and provides novel exact formulas for the recon- 
struction of the initial pressure distribution for various kinds of measurement 
setups. After the introduction of the universal back-projection algorithm intro- 
duced in [23] this goal seems superfluous, although not discussed in detail for 
sliced imaging. However, it has been shown recently by Natterer [16] that uni- 
versal back-projection is only exact for special sampling geometries. Here, for 
sliced imaging and certain sampling setups, we can indeed find mathematically 
exact universal reconstruction algorithms for arbitrary strictly convex sampling 
domains fl. 

The paper is organized as follows: In Section 2, we describe the experimental 
setup of photoacoustic sectional imaging, and we model in Section 3 various 
measurements where it is possible to derive exact reconstructions formulas for 
sectional imaging. The reconstruction formulas are then provided in Section 4. 
In the appendix, we survey some background material on the Abel transform, 
the spherical mean operator, and the Mathieu equation. 

2 Experimental Background 

Below we give an overview on photoacoustic sectional imaging, describe the 
experimental realization, and provide mathematical formulations. Opposed, we 
call conventional photoacoustic imaging with uniform illumination of the object 

non-focused. 

In general, PAI is based on the so called thermo- or photoacoustic effect. 
Laser light impinging onto a surface of an object leads to its heating and under 
the conditions of thermal and stress confinement an acoustic wave emerges from 
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the object. Mathematically, this means that the initial pressure distribution / 
is related to the absorption coefficient of the object by [3] 

= ^t^ai^Mx) (2.1) 
Op 

where x ^ . The constant prefactor consists of /3 the thermal expansivity, Cp 
the specific heat capacity and c the speed of sound. Then, /i^ is the absorption 
coefficient of the object and $ the local light fluency. Neglecting effects of light 
propagation (i.e. setting $ a constant), remains the only variable depending 
on X and containing the desired tomographic information about the object. 
Reconstruction algorithms that also take light propagation into account are 
investigated in Ref. [3] but are not focus of this work. 

It is common to classify photoacoustic measurement setups into point (see 
e.g. [11]) and integrating detector setups [2]. The focus of this paper are recon- 
struction methods for sectional imaging and various kinds of detectors. 

In conventional photoacoustic setups with point detectors, measurements 
are collected (see Fig. 1) all over a surface enclosing the object (i.e. a sphere 
or an ellipsoid), or all over a sufficiently large cylindrical surface (where the 
cross-section can be every line segment), or over a sufficiently large plane. The 
latter two detector array geometries allow for approximate reconstructions only, 
since sufficiently large means that in theory the array is considered infinitely 
large. This approximation causes the limited view problem in practical appli- 
cations. Paltauf et al. have given correction factors for some of the affected 
geometries [19]. 




Figure 1: Sketch of non- focused point detector arrays, a) closed surface, b) 
quasi-infinite plane, c) cylindrical surface. The point detectors are spread over 
the indicated surfaces (points are not drawn in the figure). The detector arrays 
b) and c) suffer from the limited view problem. 

The linear and planar detectors have to be moved tangentially to a surface 
surrounding the object (see [2]). Practically, this only allows the measurement 
devices to be aligned on a cylindrical surface (or on a plane). Experimental 
realizations of line detectors are documented for instance in [20] . 

In the following, we explain the principles of focusing detectors. The ultra- 
sonic wave is refracted by a suitable acoustic lens such that out-of-plane signals 
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Figure 2: Sketch of non-focused line detector arrays, a) quasi-infinite plane, b) 
cylindrical surface with arbitrary cross-section. The detector array a) suffers 
from the limited view problem. 

are generally weak and can be neglected. Thus only signals emerging from 
the desired imaging plane are collected at the detector. Contemporary focus- 
ing ultrasonic detectors have a spherical or cylindrical shape, thus the detector 
surface plays the role of the acoustic lens. We consider the case of a cylindri- 
cally focusing detector, which focuses into a plane. The sectional imaging can 
be improved further by illuminating the desired plane only, i.e. by cylindrical 
lenses in front of the object. Note, however, that this requires a low scattering 
coefficient of the sample, because otherwise illumination will not be restricted 
to a single plane. However, both approaches in combination provide a good 
experimental approximation to the mathematical models outlined below. The 
out-of-plane resolutions achieved are worse than the in-plane resolutions never- 
theless. For more details on focusing point detectors see [14, 21] and for focusing 
line detectors see [7, 8]. 

3 Mathematical Formulation 

We model the sectional photoacoustic imaging by assuming that the initial 
pressure distribution / : — M is perfectly focused in the illumination plane 
{x G I ^ 0}: 

/(C,^) = /(C)'5(^), CeR^zeM, (3.1) 

for some smooth function / : — )• M. The resulting pressure wave p : R'^ x 
[0,oo) — >■ R, we assume to be the solution of the linear three-dimensional wave 
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equation 



dtp{^,z;0) = 0, 
for all ^ e M^, 2; e M, and t > 0. Here, 

denotes the three-dimensional Laplacian in Euclidean coordinates. 

Our aim is to recover the function /, describing the initial pressure distri- 
bution, from certain messurements of the pressure wave p. The position of the 
detectors performing these measurements shall be given by the boundary dQ 
of a convex domain C in the illumination plane, where we additionally 
assume that / has compact support in il. 

We will consider the following four different measurement setups and derive 
reconstruction formulas. Some of the reconstruction formulas and setups have 
already been documented in the literature and are surveyed here in the general 
context. Most of the reconstruction formulas, however, are new. 

Vertical Line Detectors: The measurement data are 



mii^lt) := / p{^,z;t)dz 

J —00 

for aU ^Gdn, t>0. 



That is, in practical realizations, wc use line detectors which measure the 
overall pressure along a line orthogonal to the illumination plane. 

Point Detectors: The measurement data are 



for all ^€dn, t>0. 



That is, we use point detectors which measure the pressure on the bound- 
ary of dfl over time. This measurement geometry is used in [14, 21]. 



(3.2) 



(3.3) 



(3.4) 



For the other two measurement methods, we additionally impose that the 
domain C is strictly convex and bounded. 

Vertical Plane Detectors: The measurement data are 



mz{9\t) := I p{x\t)ds{x) 
Jp{0) 

for all 6* e S'\ t > 0, 



(3.5) 
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where P{0) C denotes the tangential plane of the cylinder dfl x R 
orthogonal to the vector {9,0), see (3.8). For a practical realization we 
use planar detectors which are moved tangentially around the object and 
measure the averaged pressure on the plane. 



Horizontal Line Detectors: The measurement data are 



7714(0; i) := / pi^,0;t)ds{^) 
JT{e) 

for all 6&S'^,t>Q, 



(3.6) 



where T{0) c R^ denotes the tangential line of dO, orthogonal to the 
vector 9, sec (3.7). This is a realization using line detectors which measure 
the overall pressure on a line tangential to dVt in the illumination plane, 
see [7, 8]. 

In those cases where the domain Q. is strictly convex and bounded, we 
parametrize the boundary dO. with the map C, : ^ dO, which associates 
to every unit vector 9 E the point (,{9) G dfl where the outward unit normal 
vector of d^l coincides with 9, see Figure 3. 



Figure 3: Definition of the point C(^)) ^ = (cos??, sini?). 

Since the tangent line T{6) of dSl at ({0) is thus by definition orthogonal 
to 9, we can define the family T(r, ^), r € M, of lines parallel to the tangent T{6) 
by 

T{r,9) = ({9) + r9 + m9^ CM.^, T{9) = T{0,9), (3.7) 

for every 9 G and r e R. 

Moreover, we introduce the family P{r,9), r e M, of planes parallel to the 
tangent plane P{9) of the cylinder 9^ x R at (C(6'), 0) by 

P{r, 6) = {T{9),0) + (0, R) C R^ P{e) = P{0, 9), (3.8) 

for every 9 G and r € R. 
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4 Reconstruction Methods 



In the following, we derive universal reconstruction formulas for photoacoustic 

sectional imaging. Conceptually, the paper is closely related to [13], where 
universal formulas for conventional photoacoustic tomography were derived in 
arbitrary geometry. However, some of the results there are implicit, requiring 
explicit knowledge of eigenfunctions, which are provided here explicitly. Even 
more, the focus of this paper is on sectional imaging, which results in different 
formulas. 

4.1 Measurements with Vertical Line Detectors 

We introduce the function 

/oo 
pitz;t)dz, ^GM.\t>0. (4.1) 
-oo 

Then the inital value problem (3.2) for the function p implies that the func- 
tion p satisfies the two-dimensional wave equation 

dttPit t) = A{p(^; t) for ah ^ G i > 

with the initial conditions 

dtp{^;0)^0 for all 

mO) = nO for ah eeM^. 

The initially three-dimensional reconstruction problem therefore reduces to the 
two-dimensional problem of calculating /(^) = p{S,',0), £. & R^, from the mea- 
surement data 

mi{^]t)=p{^;t), ^€dn,t>0. 

4.1.1 Reconstruction Formulas Based on Series Expansions 

For special domains fl, explicit reconstruction formulas are known: see the 
review [9] for V, a circle and the half-space. The derivation for the ellipse is 
published in [4]. 

• If O is the half-space e 16 > 0}, we get [10] 

/(^) = _/ / mi(A;,w)e"=«^cos(6V^j23^)dwdfc 

7^ J-oo Jk 

for every ^ e K^, where 

mi{k,uj) = - / mi{^i,0;t)e-'''^' cos{uit)dtd^i 

J-oo Jo 

is the Fourier-cosine transform of the measurement data mi . 



(4.2) 
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• If 17 = i3^(0) is the two-dimensional ball with radius R and center 0, we 
choose ip : [0,oo) x [0, 27r) — >■ R^, ip{r,(p) = (r cos ^, r sin , and get [9] 



(4.3) 



for every r G [0, oo) and (p € [0, 27r), where 



rhi{k,(jo) = — 

TT 



oo />27r 



^0 



mi{il}{R,ifi);t)e ''"^ cos{ujt) dip dt 



is the Fourier-cosine transform of the measurement data mi. Here, Jk, 
fc e No, denotes the fcth Bessel function. 



Equations (4.3) and (4.2) can be derived from formulas for the inversion of 
the spherical mean operator (see Section A. 2), and this is why these formulas 

arc typically assigned to Norton [17] and Norton & Linzer [18], although they 
considered reflectivity ultrasound imaging and here the topic is photoacoustics. 

• If O is the eUipse e j || + || < 1} with a > 6, we set 

*-:|0,»)x|0,2.)^«^ (-^MrSO- 
with the linear eccentricity e = \/o? — \P' and find [4] 



f{i^{r,v)) = —[ Y] p''/^'^\ mi(fc,a;)^fc(y;u;)rfa;, 



where ro = artanh(|) is chosen such that V(''0) € 9fl and 

mi{k,oj) = — / / mi{ip{ro,ip);t)^ki^',i^)cos{u)t)dipdt. 

Jo Jo 

Herein, the functions are for fc e Nq defined by 

^'2fc(^;w) = ccfc(^; £^) and $2fe+i(</'; w) = scfe+i(.^; S^), 

where cefe and scfc denote the Mathieu cosine and Mathieu sine functions, 
respectively, see Section A.4, and the functions Rk are the corresponding 
solutions of the radial Mathieu equation and are for k gNo given by 

R2k{r;uj) = cefe(ir; 2^) and i?2fe+i(r;w) = -isefe+i(ir; ^^). 
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4.1.2 Reduction to the Spherical Mean Operator 



Taking into account the relation (A. 8) between the solution of the two-dimen- 
sional wave equation and the spherical mean operator, the problem of deter- 
mining / from the measurements mi{^;t) = p{^;t) for ^ e di}, t > can be 
equivalently described as the problem of reconstructing / from the spherical 
mean operator A^2[/](^; of / for ^ e dfl and r e (0, oo). 

• For Cl = -B^(O) C M^, analytical reconstruction formulas have been derived 
by Finch, Haltmeier, Rakesh [6] and read as follows 





j^^ j^\{M2[mR0,r)\og\r^ - |^ - iJ^H drdsiO)^ 


and 


(4.4) 




^ {drrdrM2[f]) {R9,r)\og\r'' -\i-Re\^\drds{e). 



(4.5) 

For a general domain fi, Kunyansky reduced in [13] the reconstruction 
problem to the determination of the eigenvalues and normalized eigen- 
functions Uk, \\uk}\i = 1, of the Dirichlet Laplacian — A on f2 with zero 
boundary conditions: 

Aufe(0+Afe«fe(e) = 0, eell, (4.6) 

«fc(e) = o, ^esa (4.7) 

Indeed, if (^, 77) 1— > Ga^. (|C^??|) is a free-space rotationally invariant Green's 
function of the Helmholtz equation (4.6) and n(^) denotes the outer unit 
normal vector of dO. at ^ € dO., then 



/(0 = 27rX^M,«fe(0, 



fe=o 



(4.8) 



where 

Mfe=/ / rM2[f]{'n,r)Gx,{r){Vuk{ri),n{v)) drds{7^). 



an Jo 



00 



We thus get the initial pressure distribution / by first calculating from the 

measurements mi(^; t) = p(^; t), ^ E dfl, t > 0, with formula (A. 8) the spherical 
mean operator M2[.f ]{£,]''') of / for ^ G dVl and r € (0,oo), and then using one 
of the formulas (4.4), (4.5), or (4.8) to get /. 
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4.2 Measurements with Point Detectors 

From equation (A. 7), wc know that the solution of the initial value problem (3.2) 
can be for every x € and t > written in the form 

p{x-t)^dt[^[ f{x + y)ds{y)\. 

Parameterizing the sphere dBf{0) in cylindrical coordinates, i.e. in the form 

dBf{0) = {{Vt'^-h;^ 9,h) : h e [-t,t], 9 e S^} , 
we find for every x = {£,, z), ^ G M?, z gR, and t>0 that 

p{^, z; t) = dt (J^ y *^ j^^ /(^ + s/W^e)5{z + h)t ds{6) dhj , 

where we have used the special property of the source term (3.1). Integrating 
out the ^-distribution, we get for z G [—t, t] 

p{^, z; t) = dt (^^ y^ /(e + Vi^^e) ds{0)^ . (4.9) 

By the definition (A.4) of the spherical mean operator A42, this means 

p{i,z;t) = ^dt{M2[m;V^^^)) for zG[-t,t]. (4.10) 

From the assumption that the support of / lies completely in fl, we know that 
A^2[/](C; 0) = fiO = for ^ ^ ri. Thus, we can integrate the relation (4.10) for 
^ ^ O and find for every z G [—t, t] that 

M2[f] (C; VW^) = 2 /' p{^, z; i) di. 

J z 

Setting 2; = 0, we get for every ^ G dfl and every t > the relation 

M2[m;t) = 2 f m2{^\i)dl 
^0 

Having calculated the spherical mean of /, we can now proceed as in Sec- 
tion 4.1.2. 

4.3 Measurements with Vertical Plane Detectors 

For every 9 G S^, we define for r e M and t>0 the function 

P0{r;t) = / p{x;t)ds{x) , 
Jp(r,e) 
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where P{r,9) denotes the plane as defined in (3.8). 

Then, since the vectors {6, 0), (^-"-,0), and (0, 0, 1) form an orthonormal basis 
of and the Laplacian is rotationally invariant, we find from equation (3.2) 
that 



dttpeir;t)^ / Aa;p{C{e) + rO + ue^,z;t) dudz = drrPe{r;t) 



oo J — oo 



for every r € M and t > 0. Thus, pe solves the one-dimensional wave equation 
with the initial conditions 

dtpeir; 0) = for all r G K and 

P0{O; t) = msiO; t) for ah t > 

resulting from (3.2) and (3.5), respectively. Moreover, since / has its support 
inside il, we know that pe{r\ 0) = for r > 0. 

With d'Alembert's formula for the solution of the one-dimensional wave 
equation, we find that the unique solution for this initial value problem is given 

by 

P0{r; t) = m-i{e; -t - r) + 7713(6!; t-r), r G R, t > 0, 

where we set 7773(0; t) = for i < 0. 

Finally, we have to recover from the values ofp0, 9 G S^, the initial pressure 
distribution / from equation (3.2). We have the relation 

/•oo 

Pe{r;0)= / f{m+re + u0^)du = n[f]{r + {ae),e),0), 




where TZ denotes the Radon transform as defined in (A. 9). We can therefore 
recover / with an inverse Radon transform: 



(4.11) 



Equation (4.11) reveals an interesting property of integrating area detectors: 
For an arbitrary strictly convex measurement geometry ft, exact reconstruction 
formulas exist. This is a property which is not known for conventional and other 
photoacoustic sectional imaging technologies. 

4.4 Measurements with Horizontal Line Detectors 

For every 9 € S^,we define the function 

P0{r,z;t) = / p{^,z;t)ds{^), 

JT{r,0) 
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where T{r,d) is defined as in (3.7). Then, using that the veetors (6',0), (6*^,0), 
and (0, 0, 1) are an orthonormal basis of and that the Laplacian is rotationahy 
invariant, the initial value problem (3.2) implies that pe solves for all r,z G R 
and t > the two-dimensional wave equation 

/oo 
/:^^xP{C{0)+re + ue-^,z;t)du 
-oo 

= drrPe{r, z;t)+ d^zpeir, z; t) 

with the initial conditions 

dtP0{r,z;O) = 0, 
Pe{r,z;0) = Fe{r)6{z), Fe{r) = [ fiOdsiO, 

JT{r,0) 

for every r, 2: e R. 

Prom formula (A. 6), we see that the solution of this initial value problem 
can be written as 



27r 



for all r,z £ R and t > 0. Integrating out the 5-function, we find for every 
z e [-t, t] that 

pe{r,z;t) = —dt / —===== dp . 

Since pq is related to the measurement 7714, given by (3.6), via m^(B\i) = 
P0{O, 0; t), and since Fg{r) = for r > by the assumption that / has support 
inside O, we find with the formula (A.l) for the Abel transform in reciprocal 
coordinates that 

""("^ ') - T,"' (i' J=! ""j - i"' <^-^i«"i<«) 

where ipe{j;) — p'^Fg{~p). Switching to the reciprocal coordinate s = j and 
using the identity (A. 3), we see that this is of the form 

|™(^; -s) = -^dsisAlHi^)) = A-'[^^>e]{s) 
with ^e{^) = p-V'e(^) = Fg{—p). Thus, we can directly solve the equation for 



Fff and find 



F0{-p) = 2A[me]{^p), mg{\) = t''mi{0;t). (4.12) 
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Since we have by definition 

Fe{r) = n[f]{r + {m,()),0), 
we finally get (remembering that Fq (r) = for r > 0) 



f = 2n-^[Fl F{r,e) = 




if r < {({9), 9) , 
itr>{C{9),9). 



So, the reconstruction of / can be accomplished by an Abel transform of the 
rescaled measurements me, defined in (4.12), followed by an inverse Radon 
transform. Again, this reconstruction formula is valid for an arbitrary strictly 
convex measurement geometry fl. 

Conclusion 

In this paper we have surveyed exact reconstruction formulas for photoacoustic 
sectional imaging. All formulas arc; mathematically, analytically exact. Com- 
paring point and integrating line detectors, it is quite surprising that integrating 
area detectors allow analytical reconstructions for all strictly convex domains. 
This property has not been observed for point detectors. 



A Appendix 

A.l Abel Transform 

The Abel transform A[tp] of a smooth function tjj : 1R+ — ^ M, which decays to 
zero at oo, is defined by 

Am = r dx = 2 r -^^1= dr, y>0. 

J-oo Jy V - y 

We rewrite the Abel transform in reciprocal coordinates, so that it better fits 
in our context. Substituting y = j and r = i, we find that 

To invert the Abel transform, we remark that we have for all t; > 



J-oo Vu^ + J-oo J- 



dudx 



27r / ip'{p) dp = -2'K'4){v) 

J V 
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where we substituted x = \/ — cos((^) and u = \J — sin((/3). Therefore, 
the inverse Abel transform A~^[ip] of a function ip : IR._|_ — >■ R can be written as 



A-'my) = -^{Am'{y)- (A.2) 



Using the identity 



= dy{y^Amy)) + / , f ^ . ^'(T^^+F) dx 

J-oo ^Jx^ + t/^ 

= - V \ ^{^M+^)dx = ydy{yAmy)), 



and using (A.2) we can also write the inverse Abel transform in the form 

^-i[rV](y) = -^dyiyAmv)). (A.3) 

A.2 Circular and Spherical Means 

Let X G M" and r > 0. The spherical mean operator in M" of an integrable 
function / : M" -)> R is defined by 

Mn[f]{x;r) = [ f{x + r6)ds{6), (A.4) 

where jS'""^] denotes the area of the unit sphere S"~^ in M". 

The spherical mean value operator is closely related to the solution of the 
n-dimensional wave equation 

dttp{x;t) = Axp{x;t), 
dtp{x: 0) = 0, 
p{x;Q) = f{x) 

for all a; e M" and t > 0. More precisely, the solution p can be expressed in 
terms of the spherical mean operator of / by (see e.g. [5]) 

for all a; e M" and f > 0. In particular, we have for 

• n = 2 that the solution p of the two-dimensional wave equation can be 
calculated from the spherical means via the Abel transform (A.l): 
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where fx{^) '■= r^M2[f]{x,r), leading also to the formula 

27r \Js?(o) y/t^ - \y\^ J 

for all a; e M^, t > 0, where -Bf (0) C denotes the two-dimensional ball 
with radius t and center 0; 

• and for n = 3, we get that the solution p of the three-dimensional wave 
equation and A^sf/] are related by 

p{x;t) = dt{tM3[f]{x;t)) = dt ( / f{x + y)ds{y)] (A.7) 

for all a; e M^, f > 0, where dBf{0) denotes the boundary of the three- 
dimensional ball Bf{0) c with radius t and center 0. 

We remark that we can solve the equations (A. 5) and (A.7) for the spherical 
mean operator of /. We get for n = 2 with s = j that 

which gives us with the representation (A. 3) of the inverse Abel transform that 

M2[f]{x;r) = —A\p,]{^^) = -r ^J^^ dt, xeR^r>0, (A.8) 

where Pxij) = f^Pix'jt); and for n = 3, we find 

1 f 

M3\f](x;r) = - p(x;t)dt, xeR^,r>0. 
r Jo 

A. 3 The Radon Transform 

The Radon transform for a function £f : — R is defined by 

/oo 
g{r6 + u6^)du. (A.9) 
-oo 

The Radon transform can be inverted and we have the explicit formula 

1 f°° \ C 

T^~^[Gm = -j7r^ - drGir + {e,0,e)d9dr 
for the inverse Radon transform, see e.g. [15] . 
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A. 4 Mathieu Functions 



The Mathieu Functions are solutions of the Mathieu equation 
u"{s) + {a- 2qcos{2s))u{s) = 0. 



(A.10) 



However, we are only interested in 27r-periodic solutions. It is known, see e.g. [1] , 

that for a fixed vahic q> 0, there only exists a 27r-periodic solution of the equa- 
tion (A. 10) for a discrete set of values a S M. Conventionally, the values a 
for which an even 27r-periodic solution exists, are labeled in increasing order as 
o-niq) with the Mathieu cosine functions s i— > cen{s;q) as corresponding solu- 
tions, n e No; and the values a for which we have an odd 27r-periodic solution 
are (again in increasing order) called bn{q) with the Mathieu sine functions 
s I— > sen{s,q) as corresponding solutions, n € N. The normalization of the 
solutions ce„ and se„ is chosen to be 



Thus, since the functions are the eigenfunctions of the symmetric operator dgs + 

2gcos(2.s), the functions ^ce„(-;(j). sc„+i (•; (j), n G No, form for every 

g > a complete orthonormal system of i^([0, 2tt]). 
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